home *** CD-ROM | disk | FTP | other *** search
/ Internet Info 1994 March / Internet Info CD-ROM (Walnut Creek) (March 1994).iso / networking / applic / ntp / minimuf.c < prev    next >
C/C++ Source or Header  |  1991-09-29  |  15KB  |  437 lines
  1. /*
  2. ***********************************************************************
  3. * Program to calculate maximum usable frequency and received signal   *
  4. * strength for high-frequency radio circuits.  Uses MINIMUF 3.5.      *
  5. ***********************************************************************
  6.  
  7. input format:
  8.     first line contains five numbers:
  9.  
  10.     3 1 1 194 20
  11.  
  12.     1. output format
  13.         1 receiver power (dB/uV)
  14.         2 elevation angle (deg)
  15.         3 delay (ms)
  16.     2. month (1-12)
  17.     3. day (1-31)
  18.     4. smoothed sunspot number
  19.     5. transmitter ERP (dBW)
  20.  
  21.     second line contains transmitter coordinates and name:
  22.  
  23.     39.7000 -75.7825 W3HCF Newark 
  24.  
  25.     third and subsequent lines contain receiver coordinates
  26.     and name:
  27.  
  28.     39.6808 -75.7486 Evans Hall 
  29.     45.18  -75.45 CHU Ottawa
  30.     40.6803 -105.0408 WWV Ft Collins 
  31.     21.9906 -159.7667 WWVH Hawaii
  32.     52.3667 -1.1833 MSF Rugby
  33.  
  34. timecode transmitters, frequencies and geographic coordinates
  35.  
  36. WWV    Ft. Collins, CO    2.5, 5, 10, 15, 20 MHz    40:40N, 105:03W
  37. WWVB    Ft. Collins, CO 60 kHz                  40:40:49N, 105:02:27W
  38. WWVH    Kauai, HI       2.5, 5, 10, 15, 20 MHz  21:59:26N, 159:46:00W
  39. CHU     Ottawa, Canada  3330, 7335, 14670       45:18N, 75:45N
  40. MSF     Rugby, UK       60 kHz, 2.5, 5, 10 MHz  52:22N, 1:11W
  41.  
  42. */
  43. #include <stdio.h>
  44. #include <ctype.h>
  45. #include <math.h>
  46. /*
  47. Parameters
  48. */
  49. #define R 6371.2                /* radius of the Earth (km) */
  50. #define hE 110.                 /* height of E layer (km) */
  51. #define hF 320.                 /* height of F layer (km) */
  52. #define GAMMA 1.42              /* geomagnetic constant */
  53. #define PI 3.141592653589       /* the real thing */
  54. #define PIH PI/2.               /* the real thing/2 */
  55. #define PID 2.*PI               /* the real thing*2 */
  56. #define VOFL 2.9979250e8        /* velocity of light */
  57. #define D2R PI/180.             /* degrees to radians */
  58. #define R2D 180./PI             /* radians to degrees */
  59. #define MINBETA 5.*D2R          /* min elevation angle (rad) */
  60. #define RINZ 50.                /* receiver input impedance (ohms) */
  61. #define BOLTZ 1.380622e-23      /* Boltzmann's constant */
  62. #define NTEMP 273.              /* receiver noise temperature (deg K) */
  63. #define DELTAF 2400.            /* communication bandwidth (Hz) */
  64. #define MPATH 3.                /* multipath threshold (dB) */
  65. #define GLOSS 3.                /* ground-reflection loss (dB) */
  66. #define FMAX 8                  /* max frequencies */
  67. #define HMAX 10                 /* max hops */
  68. /*
  69. Function declarations
  70. */
  71. extern FILE *fopen();
  72. static double minimuf(void), sgn(double); void edit(int);
  73. /*
  74. Global declarations
  75. */
  76. int month;                      /* month of year (1-12) */
  77. int day;                        /* day of month*/
  78. int hour;                       /* hour of day (UTC) */
  79. double freq[FMAX];              /* frequencies (MHz) */
  80. double gain[46][FMAX];          /* antenna gain (main lobe) (dB) */
  81. double dB1;                     /* transmitter output power (dBW) */
  82. double ssn;                     /* smoothed sunspot number */
  83. double lat1;                    /* transmitter latitude (deg N) */
  84. double lon1;                    /* transmitter longitude (deg W) */
  85. char site1[30];                 /* transmitter site name */
  86. double lat2;                    /* receiver latitude (deg N) */
  87. double lon2;                    /* receiver longitude (deg W) */
  88. char site2[30];                 /* receiver site name */
  89.  
  90. double b1;                      /* transmitter bearing (rad) */
  91. double b2;                      /* receiver bearing (rad) */
  92. double d;                       /* great-circle distance (rad) */
  93. double theta;                   /* hour angle (rad) */
  94. int hop;                        /* number of ray hops */
  95. double phi;                     /* angle of incidence (rad) */
  96. FILE *fp_in;                    /* file handle */
  97. char fn_in[12] = "gain.dat";    /* antenna file name */
  98. char fq_in[12] = "ntp.dat";     /* qth file name */
  99. int flag;                       /* output format */
  100. /*
  101. Main program
  102. */
  103. main() {
  104. /*
  105. Hour variables
  106. */
  107.     int offset;             /* offset for local time (hours) */
  108.     int time;        /* local time at receiver */
  109.     char daynight[HMAX];    /* day/night indicator (jnx) */
  110.     double mufE[HMAX];      /* max E-layer muf (MHz) */
  111.     double mufF[HMAX];      /* min F-layer muf (MHz) */
  112.     double absorp[HMAX];    /* ionospheric absorption coefficient */
  113.     double dB2[HMAX];       /* receiver input power (dBuV) */
  114.     double path[HMAX];      /* path length (km) */
  115.     double beta[HMAX];      /* elevation angle (rad) */
  116.     double ant[HMAX];       /* antenna gain (dB) */
  117. /*
  118. Path variables
  119. */
  120.     double delay;           /* path distance (ms) */
  121.     double dhop;            /* hop angle/2 (rad) */
  122.     int daytime, nightime;    /* path indicators */
  123.     double psi;        /* sun zenith angle (rad) */
  124.     double xr, yr;        /* reflection zone coordinates (rad) */
  125.     double xs, ys, yp;    /* subsolar coordinates (rad) */
  126.     double fcE;        /* E-layer critical frequency (MHz) */
  127.     double fcF;        /* F-layer critical frequency (MHz) */
  128.     char cutoff;        /* layer cutoff indicator (EFM) */
  129.  
  130.     double beta1, phi1, level, muf, dist, Z, floor; /* double temporaries */
  131.     int i,j,h,n;            /* int temporaries */
  132.     float fp1, fp2;         /* float temporaries */
  133. /*
  134. Establish initial conditions
  135. */
  136.     if ((fp_in = fopen (fn_in, "r")) == NULL) {
  137.         printf ("Antenna file %s not found\n", fn_in); exit (1);
  138.         }
  139.     for (j = 0; j < FMAX; j++) {
  140.         fscanf(fp_in, " %f ", &fp1);
  141.         freq[j] = fp1;
  142.         }
  143.     for (i = 0; i < 46; i++) {
  144.         for (j = 0; j < FMAX; j++) {
  145.             fscanf(fp_in, " %f ", &fp1);
  146.             gain[i][j] = fp1;
  147.             }
  148.         }
  149.     if ((fp_in = fopen (fq_in, "r")) == NULL) {
  150.         printf ("QTH file %s not found\n", fq_in); exit (2);
  151.         }
  152.     if (fscanf(fp_in, " %i %i %i %f %f", &flag, &month, &day, &fp1, &fp2)
  153.         != 5) exit (0);
  154.     ssn = fp1; dB1 = fp2;
  155.     if (fscanf(fp_in, " %f %f %[^\n]", &fp1, &fp2, site1) != 3) exit (0);
  156.     lat1 = fp1*D2R; lon1 = -fp2*D2R;
  157.     printf("Transmitter %s\n", site1);
  158. L1:     if (fscanf(fp_in, " %f %f %[^\n]", &fp1, &fp2, site2) != 3) exit (0);
  159.     lat2 = fp1*D2R; lon2 = -fp2*D2R;
  160.     printf("Receiver %s\n", site2);
  161. /*
  162. Compute great-circle bearings, great-circle distance, min hops and path delay
  163. */
  164.     theta = lon1-lon2;
  165.     if (theta >= PI) theta = theta-PID;
  166.     if (theta <= -PI) theta = theta+PID;
  167.     d = acos(sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(theta));
  168.     if (d < 0.) d = PI+d;
  169.     b1 = acos((sin(lat2)-sin(lat1)*cos(d))/(cos(lat1)*sin(d)));
  170.     if (b1 < 0.) b1 = PI+b1;
  171.     if (theta < 0) b1 = PID-b1;
  172.     b2 = acos((sin(lat1)-sin(lat2)*cos(d))/(cos(lat2)*sin(d)));
  173.     if (b2 < 0.) b2 = PI+b2;
  174.     if (theta >= 0.) b2 = PID-b2;
  175.     hop = d/(2.*acos(R/(R+hF)));
  176.     beta1 = 0.;
  177.     while (beta1 < MINBETA) {
  178.         hop = hop+1;
  179.         dhop = d/(hop*2.);
  180.         beta1 = atan((cos(dhop)-R/(R+hF))/sin(dhop));
  181.         }
  182.     phi = R*cos(beta1)/(R+hE);
  183.     phi = atan(phi/sqrt(1.-phi*phi));
  184.     delay = 2.*hop*sin(dhop)*(R+hF)/cos(beta1)/VOFL*1e6;
  185.     printf("\nSmoothed sunspot number:%4.0f    Month:%3i     Day:%3i\n",
  186.         ssn, month, day);
  187.     printf("Power:%3.0f dBW    Distance:%6.0f km    Delay:%5.1f ms\n",
  188.         dB1, d*R, delay);
  189.     printf("Location                        Lat      Long    Azim\n");
  190.     printf("%-27s %7.2fN  %7.2fW    %3.0f\n", site1, lat1*R2D, lon1*R2D, b1*R2D);
  191.     printf("%-27s %7.2fN  %7.2fW    %3.0f\n", site2, lat2*R2D, lon2*R2D, b2*R2D);
  192.     printf("UT LT  MUF%8.1f%8.1f%8.1f%8.1f%8.1f%8.1f%8.1f%8.1f\n"
  193.         , freq[0], freq[1], freq[2], freq[3], freq[4], freq[5],
  194.         freq[6], freq[7]);
  195. /*
  196. Hour loop: This loop determines the min-hop path and next two higher-hop
  197. paths. It selects the most likely path for each frequency and calculates the
  198. receiver power. The minimum F-layer MUF is computed directly from MINIMUF 3.5
  199. and regardless of the actual number of hops, geographic coordinates or time
  200. of day.
  201. */
  202.     offset = lon2*12./PI+.5;
  203.     for (hour = 0; hour < 24; hour++) {
  204.         time = hour-offset;
  205.         if (time < 0) time = time+24;
  206.         if (time >= 24) time = time-24;
  207.         muf = minimuf();
  208.         fcF = muf * cos(phi);
  209.         printf("%2i %2i", hour, time);
  210. /* subsolar coordinates */
  211.         xs = (month-1.)*365.25/12.+day-80.;
  212.         xs = 23.5*D2R*sin(xs/365.25*PID);
  213.         ys = (hour*15.-180.)*D2R;
  214. /*
  215. Path loop: This loop determines the geometry of the min-hop path and the next
  216. two hiher-hop paths. It calculates the minimum F-layer MUF, maximum E-layer
  217. MUF and ionospheric absorption factor for each path.
  218. */
  219.         for (h = hop; h < hop+2; h++) {
  220.             daytime = 0.;
  221.             nightime = 0.;
  222.             mufE[h] = 0.;
  223.             absorp[h] = 0.;
  224.             daynight[h] = ' ';
  225. /* path geometry, min F-layer MUF */
  226.             dhop = d/(h*2.);
  227.             beta1 = atan((cos(dhop)-R/(R+hF))/sin(dhop));
  228.             phi1 = R*cos(beta1)/(R+hE);
  229.             phi1 = atan(phi1/sqrt(1.-phi1*phi1));
  230.             mufF[h] = fcF/cos(phi1);
  231. /*
  232. Hop loop: This loop calculates the reflection zones and subsolar zenith angle
  233. for each hop along the path. It then computes the minimum E-layer MUF and
  234. ionospheric absorption coeeficient for the total path.
  235. */
  236.             for (dist = dhop; dist < d; dist = dist+dhop*2) {
  237. /* reflection zone coordinates */
  238.                 xr = acos(cos(dist)*sin(lat1)+sin(dist)*cos(lat1)*cos(b1));
  239.                 if (xr < 0.) xr = PI+xr;
  240.                 xr = PIH-xr;
  241.                 yr = acos((cos(dist)-sin(xr)*sin(lat1))/(cos(xr)*cos(lat1)));
  242.                 if (yr < 0.) yr = PI+yr;
  243.                 if (theta < 0.) yr = -yr;
  244.                 yr = lon1-yr;
  245.                 if (yr >= PI) yr = yr-PID;
  246.                 if (yr <= -PI) yr = yr+PID;
  247.                 yp = ys-yr;
  248.                 if (yp > PI) yp = PID-yp;
  249.                 if (yp < -PI) yp = -PID+yp;
  250. /* E-layer MUF */
  251.                 psi = acos(sin(xr)*sin(xs)+cos(xr)*cos(xs)*cos(yp));
  252.                 if (psi < 0.) psi = PI+psi;
  253.                 Z = cos(psi);
  254.                 if (Z > 0.) {
  255.                     Z = .9*pow((180.+1.44*ssn)*Z,.25);
  256.                     if (Z < .005*ssn) Z = .005*ssn;
  257.                     }
  258.                 else Z = .005*ssn;
  259.                 Z = Z/cos(phi1);
  260.                 if (Z > mufE[h]) mufE[h] = Z;
  261. /* ionospheric absorption coeeficient */
  262.                 Z = psi;
  263.                 if (Z > 100.8*D2R) {
  264.                     Z = 100.8*D2R;
  265.                     nightime++;
  266.                     }
  267.                 else daytime++;
  268.                 Z = cos(90./100.8*Z);
  269.                 if (Z < 0.) Z = 0;
  270.                 Z = (1.+.0037*ssn)*pow(Z,1.3);
  271.                 if (Z < .1) Z = .1;
  272.                 absorp[h] = absorp[h]+Z;
  273.                 }
  274. /* path flags */
  275.             if (daytime > 0.) daynight[h] = 'j';
  276.             if (nightime > 0.) daynight[h] = 'n';
  277.             if ((daytime > 0.) && (nightime > 0.)) daynight[h] = 'x';
  278.             }
  279. /*
  280. Frequency loop: This loop calculates the receiver power for each frequency
  281. and path.  It discards paths above the minimum F-layer MUF or below the
  282. maximum E-layer MUF and selects the path with maximum power. It then
  283. calculates the combined power of the remaining paths and determines if
  284. multipath conditions exist.
  285. */
  286.         printf("%5.1f ", mufF[hop]);
  287.         for (i = 0; i < FMAX; i++) {
  288.             level = -1e6; j = 0;
  289. /* receiver power for each path */
  290.             for (h = hop; h < hop+2; h++) {
  291.                 dhop = d/(h*2.);
  292.                 switch (daynight[h]) {
  293.                     case 'n': Z = 250.; break;
  294.                     case 'j': Z = 350.; break;
  295.                     default: Z = hF;
  296.                     }
  297.                 beta[h] = atan((cos(dhop)-R/(R+Z))/sin(dhop));
  298.                 phi1 = R*cos(beta[h])/(R+hE);
  299.                 phi1 = atan(phi1/sqrt(1.-phi1*phi1));
  300.                 path[h] = 2.*h*sin(dhop)*(R+Z)/cos(beta[h]);
  301.                 dB2[h] = -1e6;
  302.                 if ((freq[i] < mufF[h]) && (freq[i] > mufE[h])) {
  303.                     Z = dB1+137.;
  304.                     Z = Z-32.44-20.*log10(path[h])-h*GLOSS;
  305.                     Z = Z-8.686*log(freq[i]);
  306.                     Z = Z-677.2*absorp[h]/cos(phi1)
  307.                         /(pow((freq[i]+GAMMA),
  308.                         1.98)+10.2);
  309.                     muf = modf(beta[h]*R2D/2., &dist);
  310.                     n = dist;
  311.                     ant[h] = gain[n][i]+muf*(gain[n+1][i]-gain[n][i]);
  312.                     Z = Z+ant[h]; dB2[h] = Z;
  313.                     if (dB2[h] > level) {
  314.                         level = dB2[h]; j = h;;
  315.                         }
  316.                     }
  317.                 }
  318. /* select path and determine multipath */
  319.             floor = 20.*log10(sqrt(4.*RINZ*BOLTZ*NTEMP*DELTAF))+120.;
  320.             muf = 0.;
  321.             if (j > 0) {
  322.                 for (h = hop; h < hop+2; h++)
  323.                     muf = muf+exp(dB2[h]/10.);
  324.                 muf = 10.*log10(muf); cutoff = ' ';
  325.                 if (dB2[j] < muf+MPATH) cutoff = 'm';
  326.                 if (dB2[j] < floor) cutoff = 's';
  327.                 switch (flag) {
  328.                     case 1: printf("%5.0f%c%1i%c",
  329.                         dB2[j],
  330.                         daynight[j], j, cutoff);
  331.                         break;
  332.                     case 2: printf("%5.1f%c%1i%c",
  333.                         beta[j]*R2D,
  334.                         daynight[j], j, cutoff);
  335.                         break;
  336.                     case 3: printf("%5.1f%c%1i%c",
  337.                         path[j]/300,
  338.                         daynight[j], j, cutoff);
  339.                     }
  340.                 }
  341.             else {
  342.                 printf("        ");
  343.                 }
  344.             }
  345.         printf("\n");
  346.         }
  347.     goto L1;
  348.     }
  349. /*
  350. MINIMUF 3.5 (From QST December 1982)
  351.  
  352. Inputs
  353. lat1 = transmitter latitude, lon1 = transmitter longitude
  354. lat2 = receiver latitude, lon2 = receiver longitude
  355. month = month of year
  356. day = day of month
  357. hour = UTC hour, ssn = sunspot number
  358.  
  359. Outputs
  360. mufF = F-layer MUF
  361. */
  362. double minimuf() {
  363.  
  364.     double MUF;
  365.     double A, B, C, D, P, Q;
  366.     double Y1, Y2;
  367.     double T, T4, T6, T9;
  368.     double G0, G1, G2, G7, G8, G9;
  369.     double K, K1, K5, K6, K7, K8, K9;
  370.     double M9, C0, U, U1, W0, L0;
  371.  
  372.     K7 = sin(lat1)*sin(lat2)+cos(lat1)*cos(lat2)*cos(lon2-lon1);
  373.     if (K7 < -1.) K7 = -1.; if (K7 > 1.) K7 = 1.; G1 = acos(K7);
  374.     K6 = 1.59*G1; if (K6 < 1.) K6 = 1.; K5 = 1./K6; MUF = 100.;
  375.     for (K1 = 1./(2.*K6); K1 <= 1.-1./(2.*K6); K1 = K1+fabs(.9999-1./K6)) {
  376.         if (K5 != 1.) K5 = .5;
  377.         P = sin(lat2); Q = cos(lat2);
  378.         A = (sin(lat1)-P*cos(G1))/(Q*sin(G1));
  379.         B = G1*K1; C = P*cos(B)+Q*sin(B)*A;
  380.         D = (cos(B)-C*P)/(Q*sqrt(1.-C*C));
  381.         if (D < -1.) D = -1.; if (D > 1.) D = 1.; D = acos(D);
  382.         W0 = lon2+sgn(sin(lon1-lon2))*D;
  383.         if (W0 < 0) W0 = W0+PID; if (W0 >= PID) W0 = W0-PID;
  384.         if (C < -1.) C = -1.; if (C > 1.) C = 1.; L0 = PIH-acos(C);
  385.         Y1 = .0172*(10.+(month-1.)*30.4+day);
  386.         Y2 = .409*cos(Y1);
  387.         K8 = 3.82*W0+12.+.13*(sin(Y1)+1.2*sin(2.*Y1));
  388.         K8 = K8-12.*(1.+sgn(K8-24.))*sgn(fabs(K8-24.));
  389.         if (cos(L0+Y2) > -.26) goto s1520;
  390.         K9 = 0.; G0 = 0.; M9 = 2.5*G1*K5; if (M9 > PIH) M9 = PIH;
  391.         M9 = sin(M9); M9 = 1.+2.5*M9*sqrt(M9);
  392.         goto s1770;
  393.  
  394. s1520:          K9 = (-.26+sin(Y2)*sin(L0))/(cos(Y2)*cos(L0)+.001);
  395.         K9 = 12.-atan(K9/sqrt(fabs(1.-K9*K9)))*7.639437;
  396.         T = K8-K9/2.+12.*(1.-sgn(K8-K9/2.))*sgn(fabs(K8-K9/2.));
  397.         T4 = K8+K9/2.-12.*(1.+sgn(K8+K9/2.-24.))*sgn(fabs(K8+K9/2-24.));
  398.         C0 = fabs(cos(L0+Y2));
  399.         T9 = 9.7*pow(C0,9.6); if (T9 < .1) T9 = .1;
  400.         M9 = 2.5*G1*K5; if (M9 > PIH) M9 = PIH;
  401.         M9 = sin(M9); M9 = 1.+2.5*M9*sqrt(M9);
  402.         if (T4 < T) goto s1680;
  403.         if ((hour-T)*(T4-hour) > 0.) goto s1690;
  404.         goto s1820;
  405.  
  406. s1680:          if ((hour-T4)*(T-hour) > 0.) goto s1820;
  407. s1690:          T6 = hour+12.*(1.+sgn(T-hour))*sgn(fabs(T-hour));
  408.         G9 = PI*(T6-T)/K9; G8 = PI*T9/K9; U = (T-T6)/T9;
  409.         G0 = C0*(sin(G9)+G8*(exp(U)-cos(G9)))/(1.+G8*G8);
  410.         G7 = C0*(G8*(exp(-K9/T9)+1.))*exp((K9-24)/2.)/(1.+G8*G8);
  411.         if (G0 < G7) G0 = G7;
  412.         goto s1770;
  413.  
  414. s1770:          G2 = (1.+ssn/250.)*M9*sqrt(6.+58.*sqrt(G0));
  415.         G2 = G2*(1.-.1*exp((K9-24.)/3.));
  416.         G2 = G2*(1.+(1.-sgn(lat1)*sgn(lat2))*.1);
  417.         G2 = G2*(1.-.1*(1.+sgn(fabs(sin(L0))-cos(L0))));
  418.         goto s1880;
  419.  
  420. s1820:          T6 = hour+12.*(1.+sgn(T4-hour))*sgn(fabs(T4-hour));
  421.         G8 = PI*T9/K9; U = (T4-T6)/2.; U1 = -K9/T9;
  422.         G0 = C0*(G8*(exp(U1)+1.))*exp(U)/(1.+G8*G8);
  423.         goto s1770;
  424.  
  425. s1880:          if (MUF > G2) MUF = G2;
  426.         }
  427.     return (MUF);
  428.     }
  429. /*
  430. sgn function (like BASIC)
  431. */
  432. double sgn(double x) {
  433.     if (x > 0.) return (1.);
  434.     if (x < 0.) return (-1.);
  435.     return (0);
  436.     }
  437.